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Abstract 

We consider a fourth order extension to MacCor- 
mack’s scheme. The original extension was fourth 
order only for the inviscid terms but was second or- 
der for the viscous terms. We show how to modify 
the viscous terms so that the scheme is uniformly 
fourth order in the spatial derivatives. Applica- 
tions are given to some boundary layer flows. In 
addition, for applications to shear flows the effect 
of the outflow boundary conditions are very impor- 
tant. We compare the accuracy of several of these 
different boundary conditions for both boundary 
layer and shear flows. Stretching at the outflow 
usually increases the oscillations in the numerical 
solution but the addition of a filtered sponge layer 
(with or without stretching) reduces such oscilla- 
tions. The oscillations are generated by insufficient 
resolution of the shear layer. When the shear layer 
is sufficiently resolved then oscillations are not gen- 
erated and there is less of a need for a nonreflecting 
boundary condition. 


1 Basic Scheme 

We consider the equation 

Ut ~b fx = (i>(x)tl x ) x (1) 

The original scheme suggested by Gottlieb and 
Turkel [5] was 
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u(1) = u " + ^x { ~ F ” +2 + SF ' n+l ~ 7Fn 


where 


F? = [-/. + b t 


(Ui - tti-l). 


Ax 


( 2 ) 


with a second stage 

n + l _ < + + _^L( 7 jf> - 8i^ + F^\) 

2 12 Ax' ' * 1 ' 2 


if 1 = i-/.+t» ] (1) 

This scheme is second order in time and third 
order in space if b(x) = 0, but is only second order 
in space for the viscous term. We label the above 
scheme FB since it has a forward difference on the 
first stage and a backward difference in the second 
stage. We next present the BF variant as 

tit 1 ) = „» + £-{F?_ 3 - 8 jFJLj + 7 F?) 
bAx 

with a second stage 

«?« - - 8F -t + 

where 

It is shown in [5] that by rotating these two 
variants (i.e. alternating the order of the sweeps) 
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one can obtain fourth order accuracy in space when 
b(x) = 0. However, one still obtains only second 
order accuracy for the viscous terms. When the 
Reynolds number is large then the dominant error 
comes from the inviscid portion and the total error 
is essentially fourth order in space. However, there 
are cases when one wishes the spatial error to be 
uniformly fourth order. We propose to accomplish 
this by a small modification to (2). The forward 
step is modified to : 

F? = [-/, + i(-6< + 8 

A similar modification is used on the backward 
differences of F. The stencil is now widened in 
that point i in the predictor requires a centered 
five point stencil as does the corrector. This wider 
stencil occurs only because of the coefficient b. If 
b is a function only of x then this extended stencil 
is not of any importance. For the Navier-Stokes 
equations the function b is essentially the viscosity 
ft. If one uses Sutherland’s law to compute then 
b is a nonlinear function and so one needs to store 
data in artificial cells for both the predictor and 
corrector at all boundaries. Thus, the way that 
the scheme is now used is: update the solution at 
all interior points. Use third order extrapolation 
to define the fluxes at two artificial points outside 
all the boundaries (see [1]). Apply the corrector 
at all interior points and then again extrapolate to 
the artificial points. At all solid boundaries the 
boundary conditions are imposed before the ex- 
trapolation. We stress that the extrapolation of 
the fluxes is identical to using one sided differences 
at the next update. The use of extrapolation rather 
than one sided differences is done only to improve 
the vectorization of the algorithm. A Taylor series 
expansion verifies that after rotating between the 
two variants the scheme is now fourth order accu- 
rate for both the hyperbolic and parabolic portions 
of the scheme. 

2 Boundary Conditions 

One of the main difficulties in solving the Navier- 
Stokes equations is the treatment of the outflow 
boundary condition. We are mainly interested in 
the viscous but high Reynolds number flow. The 
first question is whether the boundary treatment 
should be based on the Euler equations or the Navier- 
Stokes equations. The difference between the two 
approaches is not just the type of boundary con- 
ditions but even the number of boundary condi- 
tions that needs to be externally given. For invis- 
cid flow, when the flow is subsonic one boundary 


condition needs to be specified corresponding to 
an incoming acoustic wave. For supersonic flow 
no boundary conditions are specified. For viscous 
the flow the equations are no longer hyperbolic but 
rather incompletely parabolic. For subsonic flow 
four conditions need to be specified (in two dimen- 
sions), while for supersonic flow three conditions 
need to be specified. Moreover, for the inviscid 
(hyperbolic) case the needed boundary conditions 
should be of Dirichlet type, i.e. a combination of 
the dependent variables is specified. For the vis- 
cous problem (parabolic) the specified boundary 
condition can be either of Dirichlet type or Neu- 
mann type (combinations of normal derivatives) or 
a combination of both of these (Robin type). Many 
codes use inviscid type boundary conditions. This 
is based on the assumption that the flow in the 
far field is essentially inviscid because of the high 
Reynolds number and the lack of physical bound- 
ary layers. 

To be more specific we need to consider dif- 
ferent types of configurations. For boundary layer 
flows one needs to distinguish between the portion 
of the computational domain inside and outside the 
boundary layer. Outside the boundary layer one 
may be able to use inviscid type boundary con- 
ditions. Inside the boundary layer the pressure 
should be specified. Gustaffson [7] has shown that 
the problem is also well posed if one extrapolates 
all the variables inside the boundary layer. For 
external flow about an airfoil some codes use in- 
viscid type conditions while others extrapolate all 
the variables. Though the flow seems to be inviscid 
in the far field nevertheless viscous effects persist 
in the wake region. Thus, for example, there is a 
velocity defect no matter how far one goes down- 
stream and the integral of this is constant. Thus, 
as one proceeds further downstream the defect lo- 
cally gets smaller but is spread over a larger region. 
This should affect the appropriate boundary condi- 
tions. In this paper we will also consider free shear 
flows. Here too viscous effects should be impor- 
tant near the shear layer even far downstream. A 
further complication that is most pronounced for 
shear flows is that one does not know the solution 
downstream and therefore one cannot impose any 
type of Dirichlet boundary condition. Frequently, 
there is a significant spreading of the shear layer 
and so one does not know in advance even where 
the shear layer will intersect the outflow boundary. 
Furthermore, many theories expand the solution 
about a constant pressure in the far field and so 
obtain a differential equation for the pressure de- 
viation. For shear flows the pressure differs on the 
two sides of the shear and so the pressure is not con- 
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stant in the fai field except if the far field boundary 
is extremely far away which is not computationally 
practical. In other words, some of the nonreflect- 
ing boundary conditions that have been proposed 
in the literature are based on suppositions of the 
form of the outgoing wave, e.g. a plane wave or 
a spherical wave. However, these assumptions are 
not valid for shear flows. 

In spite of all these dangers we shall consider 
characteristic-like boundary conditions at outflow 
and so the number of boundary conditions will be 
given by the inviscid theory. Nevertheless we shall 
see that viscous effects are at least partially ac- 
counted for. 

1. BC1: The simplest approach is to freeze 

the characteristic variables normal to the sur- 
face and to specify the incoming characteris- 
tic variable and to extrapolate the outgoing 
variables. For a one dimensional hyperbolic 
system one can show that such a procedure 
is well posed [6]. This approach was used by 
many authors in the past. 

2. BC2: An improved version of this system is 
to use differential equations that correspond 
to these characteristic variables. Thus, for 
the acoustic waves one needs differential equa- 
tions for p t ± petit where u is the velocity 
component normal to the boundary. For the 
shear wave we need Vt where v is tangen- 
tial to the boundary and finally p t - c 2 p t for 
the entropy variable. Whenever the bound- 
ary condition is not specified but free to float 
then the appropriate characteristic variable is 
updated by the partial differential equation. 
In order to avoid one sided differences the 
fluxes are extrapolated outside the domain 
to artificial points. Whenever the appropri- 
ate combination is specified then we replace 
this by specifying the combination of the time 
derivatives. We can describe this as 

p t - pcu t = R\ 
p t + pcu t - R 2 
Pt ~ c 2 p t - Rz 
v t = R 4 

where Ri is determined by which variables 
are specified and which are not. Whenever, 
the combination is not specified that Ri is 
just those spatial derivatives that come from 
the Navier-Stokes equations. Thus, Ri con- 
tains viscous contributions even though the 


basic format is based on inviscid characteris- 
tic theory. In implementing these differential 
equations we convert them to conservation 
variables p, m — pu y n — pv, E. Assuming an 
ideal gas we then have 

u 2 -f- v 2 N 

p t = (7 - 1 ){E t + r Pt - um t ~ vn t ) 


rrif upt 



For subsonic flow, the immediate generaliza- 
tion of the first method is to set R\ = 0 
and to calculate R 2 , Rz , R 4 from the Navier- 
Stokes equations. 

3. BC3: We next consider improvements for 

the last boundary condition. These improve- 
ments all leave R 2 , Rz, R 4 as before from the 
Navier-Stokes equations. The changes are all 
to the incoming wave, Giles [4] (and later 
Kroner [10]) added some y (tangential) space 
derivatives to the first equation. Thus the 
equation for R\ is replaced by 

Pt - pCUt + UpCVy 4* v(p y — pCUy) — 0 

As before all the derivatives can be trans- 
formed to derivatives of conservation vari- 
ables. 

4. BC4: Based on an asymptotic expansion of 
the wave equation Bayliss and Turkel [2] de- 
rived the following nonreflecting condition to 
replace R\. 

Define d 2 ~ + y 2 where M is the Mach 

number, x, y are the physical locations of 
the boundary point relative to some source, 
usually the inflow. Then 

vt - - uv ^ 

+ c ~^ 1 {v ~ Poo) 

= 0 (3) 

As stated before one frequently does not know 
Poo and so we shall simply ignore the last 
term in this equation. 

We note that this equation does not have the 
form of the first equation for R\ y i.e. it is not 
an equation for a characteristic variable. 
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5. BC5: For m&ny cases the domain is much 
longer than it is high. In this case we can 
ignore y relative to x. Hence, we assume | ~ 
\/l - If 5 , | ~ 0 . Then (3) simplifies to 

Pt - pc(u t -uv y ) = 0 (4) 

With this simplification we get an equation 
for the characteristic variable in a form simi- 
lar to that proposed by Giles but not identical 
with his condition. 

Besides nonreflecting boundary conditions based 
on the partial differential equation there are other 
approaches to removing the difficulties associated 
with far field boundaries. One popular approach is 
to use a ‘sponge layer 5 . The idea behind this ap- 
proach is to introduce a domain between the region 
of interest and the actual outflow boundary. In this 
region one changes the differential equation in some 
manner so as dissipate the energy of the system or 
else prevent waves from traveling back into the do- 
main. Since, the equations have been changed the 
solution in this sponge layer has no physical rele- 
vance. Hopefully, the solution in the rest of the do- 
main is physically relevant and is not contaminated 
with any false reflections from the boundary or the 
sponge layer itself. In this approach the boundary 
condition at the outflow boundary is irrelevant. It 
may introduce perturbations but these are elimi- 
nated in the sponge layer. For this approach to be 
effective the sponge layer must be small relative to 
the rest of the domain. A sponge layer approach 
suggested by Colonius et. al [3] is to stretch the 
grid and then filter the solution near the outflow 
boundary. They used stretching throughout the 
domain and applied filter near the outflow (in the 
sponge layer). Karni [9] used sponge layer with 
modified governing equations to accelerate conver- 
gence to steady state. 

3 Results 

In this section we check on the improvements to 
the 2-4 scheme presented above. We consider flow 
over a flat plate with M = 0.25, Re = 100 and 
Pr — 0.72. We use a uniform cartesian mesh in 
both the directions. In figure 1 we compare the 
original method which is fourth order accurate only 
for the inviscid terms but second order for the vis- 
cous part with the improved method which is uni- 
formly fourth order accurate. In this figure we plot 
the u component of the velocity versus the nor- 
malized y distance ( 17 ). We consider two meshes. 
The finer mesh is 400 x 80 and the coarser mesh 


is 200 x 40. On the finer mesh both the original 
scheme and the improved method give similar re- 
sults. On the coarser mesh we clearly see the im- 
provement that comes from using a fourth order 
accurate treatment of the viscous terms. 

We next consider a shear flow [ 8 ] modelling a 
planar jet. All the boundary conditions used in 
this study are based on the assumption of an out- 
going plane wave expansion or else a circular wave 
expansion. Neither of these assumptions is valid 
for a jet geometry. Jet acoustics presents an addi- 
tional difficulty for far field boundaries. The stan- 
dard assumption is that the further out the artifi- 
cial boundary is, the more accurate the results are 
since it better simulates an infinite domain. How- 
ever, in many cases the flow is spatially unstable 
and as the length of the domain in the x dimen- 
sion gets larger, the waves become more unstable 
and grow until nonlinear behavior either saturates 
the wave growth or else the code stops running be- 
cause of an explosive nonlinear growth. Neverthe- 
less, we shall use these boundary conditions as the 
best available. The basic mesh is uniform in the x 
direction and stretched in the y direction. The size 
of the standard domain is 100 x 5 with 600 x 60 grid 
points. We also consider the addition of a layer be- 
yond x = 100 . This additional grid which we call 
a sponge layer is just 60 additional grid points be- 
yond x = 100 until x = 133 with a stretched mesh. 
We use a stretching proportional to x 1 5 . For ease 
of comparison all figures use the same scaling and 
show 160 radii in the axial and 5 radii in trans- 
verse direction respectively. The figures are shown 
after a nondimensional time of 500 which requires 
about a hundred thousand time steps. The solu- 
tion should be approaching a steady state. 

We now consider a case where both the inner 
and outer flows are subsonic, with the jet at M=0.8 
and the outer flow at M=0.28. Hence, the nonre- 
flecting boundary conditions are applied along the 
entire outflow line. The Reynolds number for this 
flow is 10,200 based on the (planar) jet radius. We 
compare the different boundary conditions and the 
effect of sponge layer with and without filtering. 
Figure 2 a shows the contour plot of the vorticity 
for the standard case with a finer mesh in the y di- 
rection and an extended, though uniform, domain 
in the x direction, so that the mesh is now 900 x 100 
and x < 150. We do not have a way of quantita- 
tively comparing these solutions. Hence, we shall 
content ourselves with qualitatively comparing the 
solutions with this finer grid solution. We now con- 
sider the outflow boundary conditions BC 2 - BC 5 
respectively. In figures 2 b - 2 e we plot the vortic- 
ity for the same physical case as before but with 
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the standard mesh, i.e. 600 x 60 and x < 100. 
The stretching in the y direction is stronger in the 
finer mesh (figure 2a case) than the standard mesh 
(fig. 2b-2e). All the boundary conditions on the 
same mesh give similar solutions. Note, that in 
all these cases the vorticity is visible much closer 
to the inflow than when the outer boundary was 
at x = 150. This demonstrates that for jet flows 
the position of the outer boundary is more impor- 
tant that the details of the nonreflecting boundary 
condition. Among the various boundary conditions 
the condition BC2 is slightly worse than the others. 

It is our observation that the mesh density in the 
y direction has a greater effect than the difference 
between the various boundary conditions. 

We next consider the introduction of the sponge 
layer with the nonuniform increasing grid. In this 
layer we add 60 grid points at the outflow bound- 
ary to the previous 600 points. The mesh in the 
sponge layer is stretched in the x direction using 
the stretching function x 1 5 . Stretching in y is the 
same as the previous 600 x 60 grid cases. The physi- 
cal mesh is now x < 150. Because of this stretching 
the resolution of the vortices near the outer bound- 
ary is severely reduced. Now, the effect of the dif- 
ferent boundary conditions is negligible. Hence, in 
figure 3 we plot the results obtained with the Giles 
boundary condition. The vortex growth is slightly 
delayed compared with figure 2.3 while the vortices 
near x — 100 are stronger. As expected the vortices 
near the new outer boundary are washed out. 

We next introduce the filter given below in last 
160 points near the outflow boundary, i.e. the 60 
points of the sponge layer and an additional 100 
interior points. The filter is given by 

fi = .625/* + .25(/t+i + /t-i) - .0625(/i+2 + fi- 2 ) 

This filter is 4th order accurate for a uniform mesh. 
The filter function used by Colonius et. al. [3] has 
a similar form but variable coefficients. We plot 
the vorticity field for this case in figure 4. Again, 
the numerical solution is independent of the far 
field boundary condition and so in figure 4 we only 
consider the case using the Giles boundary condi- 
tion. We see that the stretched mesh coupled with 
a larger domain delays the growth of the vorticities. 
Hence, the beginning of the vortex growth is closer 
to figure 2a than the previous cases. In the far field 
the vortex growth has been completely destroyed 
by the filtering. Hence, we only expect accurate 
solutions for x < 100. 

We finally consider a flow with the jet enter- 
ing at M— 1.5 while the outer region has an inflow 
of M=0.53. Since the inner region has a super- 


sonic flow at the exit all the variables are extrap- 
olated. The nonreflecting boundary conditions are 
used only in the outer region. The Reynolds num- 
ber based on the jet radius is about 6.37 x 10 3 . 
In figure 5a we consider the fine resolution grid of 
900 x 100 as in the subsonic case. In this case there 
is a much lower growth rate than in the subsonic 
case and the vortices are barely forming at x = 150. 
In figures 5b-5e, we compare the solutions of the 
flow computed with the various outflow boundary 
conditions BC2-5. If we reduced the mesh in the 
y direction to 60 points, as in the subsonic case, 
then the mesh is not fine enough to allow for vor- 
tex growth. Hence, for the smaller domain we shall 
consider a mesh of 600 x 100 with x < 100. In this 
mesh we use the same stretching in the y direction 
as in the subsonic cases (fig. 2b- 2e). In this case 
the Giles boundary condition seems to be slightly 
worse than the others and generates more, false, 
vorticity. On the other hand if we add a sponge 
layer with a stretched mesh, as before, the only 
reasonable solutions are given by the Giles bound- 
ary condition. Adding the filter, as given above, 
eliminated all the vortices and the solution is es- 
sentially independent of the outer boundary treat- 
ment. This is the best solution but is special for 
this case. Hence, the supersonic case is less useful 
for comparing treatments of the outer boundary. 
All the boundary conditions considered here are 
based on the inviscid case even though we compute 
the full Navier-Stokes equations. Thus, we assume 
that the viscosity is negligible at the outflow. 

4 Conclusions 

We have shown how to modify the fourth order ex- 
tension of MacCormack’s scheme so that it is uni- 
formly fourth order accurate for both the inviscid 
and viscous portions of the flux. This results in in- 
creased accuracy for boundary layer type flows at 
local unit Reynolds numbers of 1000 and lower. 

We next compared several boundary conditions 
for jet flow. In all cases it is best to use the par- 
tial differential equations at the outflow boundary 
itself. This can be accomplished by either using 
one sided differences or else some extrapolation to 
artificial points beyond the boundary. This extrap- 
olation is equivalent to a one sided difference for- 
mula within the differential operator solver. One 
then takes combinations of these updated differ- 
ences together with a radiation boundary condition 
to form the final updated solution. The combi- 
nation of given boundary conditions and equation 
solvers is determined by characteristic theory for 
the inviscid portion of the system. Nevertheless, 
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updating the complete equations introduces some 
viscous information into the procedure. 

Though the differences were not large, the best 
results were obtained with the radiation boundary 
condition of Giles or else the simplification of the 
boundary condition of Bayliss and Turkel. An al- 
ternative is to stretch the mesh at the outflow and 
then to use the obtained solution only in the orig- 
inal domain. This requires extra storage and com- 
puter time but yields somewhat better solutions. 
In this case all the boundary treatments at the out- 
flow boundary give essentially identical results as 
the major effect is the stretching of the mesh. In- 
troduction of a filter in this sponge layer with the 
expanding mesh smooths out all significant features 
of the solution in the far field. 
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Figure 2a: Subsonic jet simulation with high resolution in large domain 
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Figure 2b: Subsonic jet with BC2 boundary condition 
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Figure 2c: Subsonic jet with Giles’ boundary condition 
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Figure 2d: Subsonic jet with Bayliss-Turkel boundary condition 
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Figure 2e: Subsonic jet with simplified Bayliss-Turkel condition 
Figure 2: Comparison of boundary treatments for subsonic flows 
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Figure 3: Subsonic jet with sponge layer and without filter 


VORTICITY MAGNITUDE 



Figure 4: Subsonic jet with sponge layer and filter 
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Supersonic jet simulation with high resolution in large domain 
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Figure 5b: Supersonic jet with BC2 boundary condition 
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Figure 5c: Supersonic jet with Giles’ boundary condition 
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Figure 5d: Supersonic jet with Bayliss-Turkel boundary condition 
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Figure 5e: Supersonic jet with simplified Bayliss Turkel boundary condition 


Figure 5: Comparison of boundary treatments for supersonic flows 
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of several of these different boundary conditions for both boundary layer and shear flows. Stretching at the outflow 
usually increases the oscillations in the numerical solution but the addition of a filtered sponge layer (with or without 
stretching) reduces such oscillations. The oscillations are generated by insufficient resolution of the shear layer. When 
the shear layer is sufficiently resolved then oscillations are not generated and there is less of a need for a nonreflecting 
boundary condition. 
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